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This review explains in a self-contained way the properties of random Boolean networks and 
their attractors, with a special focus on critical networks. Using small example networks, analytical 
calculations, phenomenological arguments, and problems to solve, the basic concepts are introduced 
and important results concerning phase diagrams, numbers of relevant nodes and attractor properties 
are derived 
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1. INTRODUCTION 

Random Boolean networks (RBNs) were introduced in 
1969 by S. Kauffman 0,0 MS R simple model for gene reg- 
ulatory networks. Each gene was represented by a node 
that has two possible states, "on" (corresponding to a 
gene that is being transcribed) and "off" (corresponding 
to a gene that is not being transcribed). There are al- 
together N nodes, and each node receives input from K 
randomly chosen nodes, which represent the genes that 
control the considered gene. Furthermore, each node is 
assigned an update function that prescribes the state of 
the node in the next time step, given the state of its in- 
put nodes. This update function is chosen from the set 
of all possible update functions according to some prob- 
ability distribution. Starting from some initial configu- 
ration, the states of all nodes of the network are updated 
in parallel. Since configuration space is finite and since 
dynamics is deterministic, the system must eventually re- 
turn to a configuration that it has had before, and from 
then on it repeats the same sequence of configurations 
periodically: it is on an attractor. 

S. Kauffman focussed his interest on critical networks, 
which are at the boundary between frozen networks with 
only very short attractors and chaotic networks with at- 
tractors that may include a finite proportion of state 
space. He equated attractors with cell types. Since each 
cell contains the same DNA (i.e., the same network), cells 
can only differ by the pattern of gene activity. Based on 
results of computer simulations for the network sizes pos- 
sible at that time, S. Kauffman found that the mean num- 
ber of attractors in critical networks with K = 2 inputs 
per node increases as s/N. This finding was very satis- 
fying, since the biological data available at that time for 
various species indicated that the number of cell types is 
proportional to the square root of the number of genes. 
This would mean that the very simple model of RBNs 
with its random wiring and its random assignment of up- 
date functions displays the same scaling laws as the more 
complex reality. The concept of universality, familar from 
equilibrium critical phenomena, appeared to work also 
for this class of nonequilibrium systems. Kauffman found 
also that the mean length of attractors increases as \J~N . 

Today we know that the biological data and the com- 
puter simulation data arc both incorrect. The sequencing 



of entire genomes in recent years revealed that the num- 
ber of genes is not proportional to the mass of DNA (as 
was assumed at that time), but much smaller for higher 
organisms. The square-root law for attractor numbers 
and lengths in RBNs survived until RBNs were studied 
with much more powerful computers. Then it was found 
that for larger N the apparent square-root law does not 
hold any more, but that the increase with system size is 
faster. The numerical work was complemented by several 
beautiful analytical papers, and today we know that the 
attractor number and length of K = 2 networks increases 
with network size faster than any power law. We also 
know that, while attractor numbers do not obey power 
laws, other properties of critical RBNs do obey power 
laws. 

It is the purpose of this review to explain in an under- 
standable and self-contained way the properties of RBNs 
and their attractors, with a special focus on critical net- 
works. To this aim, this review contains examples, short 
calculations, phenomenological arguments, and problems 
to solve. Long calculations and plots of computer simu- 
lation data were not included and are not necessary for 
the understanding of the arguments. The readers will 
also benefit from consulting the review [|[, which, while 
not containing the more recent findings, covers many im- 
portant topics related to Boolean networks. 

Boolean networks are used not only to model gene reg- 
ulation networks, but also neural networks, social net- 
works, and protein interaction networks. The structure 
of all these networks is different from RBNs with their 
random wiring and random assignment of update func- 
tions, and with the same number of inputs for every node. 
Nevertheless, understanding RBNs is a first and impor- 
tant step on our way to understanding the more complex 
real networks. 



2. MODEL 

A random Boolean network is specified by its topology 
and its dynamical rules. The topology is given by the 
nodes and the links between these nodes. The links arc 
directed, i.e., they have an arrow pointing from a node 
to those nodes that it influences. The dynamical rules 
describe how the states of the nodes change with time. 
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The state of each node is "on" or "off" , and it is deter- 
mined by the state of the nodes that have links to it (i.e., 
that are its inputs). In the following, we first describe 
the topology, and then the dynamics of RBNs. 



A. Topology 

For a given number N of nodes and a given number 
K of inputs per node, a RBN is constructed by choosing 
the K inputs of each node at random among all nodes. 
If we construct a sufficiently large number of networks in 
this way, we generate an ensemble of networks. In this 
ensemble, all possible topologies occur, but their statis- 
tical weights are usually different. Let us consider the 
simplest possible example, N — 2 and K = 1, shown 
in Figure I2A1 There are 3 possible topologies. Topolo- 




(a) (b) (<0 



FIG. 1: The possible topologies for N = 2 and K = 1. 

gies (a) and (b) have each the statistical weight 1/4 in 
the ensemble, since each of the links is connected in the 
given way with probability 1/2. Topology (c) has the 
weight 1/2, since there are two possibilities for realizing 
this topology: either of the two nodes can be the one 
with the self-link. 

While the number of inputs of each node is fixed by 
the parameter K, the number of outputs (i.e. of outgo- 
ing links) varies between the nodes. The mean number of 
outputs must be K, since there must be in total the same 
number of outputs as inputs. A given node becomes the 
input of each of the N nodes with probability K/N. In 
the thermodynamic limit N — > oo the probability distri- 
bution of the number of outputs is therefore a Poisson 
distribution 

P out (k) = ^e~ K . (1) 



B. Update functions 

Next, let us specify the dynamics of the networks. 
Each node can be in the state crj = 1 ("on") or in the 
state <7i — ("off"), where i is the index of the node. 
The N nodes of the network can therefore together as- 
sume 2 N different states. An update function specifies 
the state of a node in the next time step, given the state 
of its K inputs at the present time step. Since each of the 
K inputs of a node can be on or off, there are M = 2 K 
possible input states. The update function has to specify 



TABLE I: The 4 update functions for nodes with 1 input. The 
first column lists the 2 possible states of the input, the other 
columns represent one update function each, falling into two 
classes. 
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TABLE II: The 16 update functions for nodes with 2 inputs. 
The first column lists the 4 possible states of the two inputs, 
the other columns represent one update function each, falling 
into four classes. 
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the new state of a node for each of these input states. 
Consequently, there are 2 M different update functions. 

Table U lists the 4 possible update functions for K = 
1. The first two functions are constant, or "frozen", i.e. 
the state of the node is independent of its inputs. The 
other two functions change whenever an input changes, 
i.e., they are reversible. The third function is the "copy" 
function, the fourth is the "invert" function. 

Table |TT] lists the 16 possible update functions for 
K = 2. There are again two constant and two reversible 
functions. Furthermore, there are canalizing functions. 
A function is canalyzing if at least for one value of one of 
its inputs the output is fixed, irrespective of the values of 
the other inputs. The first class of canalyzing functions 
do not depend at all on one of the two inputs. They 
simply copy or invert the value of one of the inputs. In 
Table [TTl these are the C\ functions. The second class of 
canalyzing functions has three times a 1 or three times 
a in its output (the C 2 functions). For each of the two 
inputs there exists one value that fixes the output irre- 
spective of the other input. In fact, constant functions 
can also be considered as canalyzing functions, because 
the output is fixed for any value of the inputs. 

Each node in the network is assigned an update func- 
tion by randomly choosing the function from all possible 
functions with K inputs according to some probability 
distribution. The simplest probability distribution is a 
constant one. For K = 2 networks, each function is then 
chosen with probability 1/16. In the previous section, we 
have introduced the concept an ensemble of networks. If 
we are only interested in topology, an ensemble is de- 
fined by the values of N and K. When we want to study 
network dynamics, we have to assign update functions 
to each network, and the ensemble needs to be specified 
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by also indicating which probability distribution of the 
update functions shall be used. If all 4 update functions 
are allowed, there are 36 different networks in the en- 
semble shown in Figure [2A"1 For topologies (a) and (b), 
there are 10 different possiblities to assign update func- 
tions, for topology (c) there are 16 different possibilities. 
The determination of the statistical weight of each of the 
36 networks for the case that every update function is 
chosen with the same probability is left to the reader.... 

In the following we list several frequently used proba- 
bility distributions for the update functions. Throughout 
this article, we will refer to these different "update rules" . 

1. Biased functions: A function with n times the out- 
put value 1 and M — n times the output value is 
assigned a probability p n (l — p) AI ~ n . Then the two 
frozen functions in table HI] have the probabilities 
p 4 and (1 — p) 4 , each of the C\ functions and of the 
reversible functions has the probability p 2 (l — p) 2 , 
and the C2 functions have the probabilities p(l— p) 3 
and p 3 ( 1 — p) . For the special case p = 1 /2 , all func- 
tions have the same probability 1/16. 

2. Weighted classes: All functions in the same class 
are assigned the same probability. K = 1 net- 
works are most interesting if the two reversible 
functions occur with probability 1/2 each, and the 
two constant functions with probability 0. In gen- 
eral K = 1 networks, we denote the weight of the 
constant functions with 8. An ensemble of K = 2 
networks is specified by the four parameters a, [3, 
7, and S for the weight of Ci, reversible, C2 and 
frozen functions. The sum of the four weights must 
be 1, i.e., l = a + /3 + 7 + < ! >. 

3. Only canalyzing functions are chosen, often includ- 
ing the constant functions. This is motivated by 
the finding that gene regulation networks appear 
to have many canalyzing functions and by consid- 
erations that canalyzing functions are biologically 
meaningful 0, [|| . Several authors create canalyz- 
ing networks using three parameters Q . One input 
of the node is chosen at random to be a canalyzing 
input. The first parameter, 77, is the probability 
that this input is canalyzing if its value is 1. The 
second parameter, r, is the probability that the out- 
put is 1 if the input is on its canalyzing value. The 
third parameter, p, assigns update functions for the 
K — 1 other inputs according to rule 1 (biased func- 
tions), for the case that the canalyzing input is not 
on its canalyzing value. (This notation is not uni- 
form throughout literature. For instance, in Q, the 
second and third parameter are named pi and p2-) 

4. Only threshold functions are chosen, i.e. the up- 
date rule is 

. + 11 = /lifEf =1 M2.,-l) +/l )>0 



The couplings are zero if node i receives no input 
from node j, and they are ±1 with equal probabil- 
ity if node j is an input to node i. Negative cou- 
plings are inhibitory, positive couplings are excita- 
tory. The parameter h is the threshold. Threshold 
networks are inspired by neural networks, but they 
are also used in some models for gene regulation 
networks 

5. All nodes are assigned the same function. The 
network is then a cellular automaton with random 
wiring. 



C. Dynamics 

Throughout this paper, we only consider the case of 
parallel update. All nodes are updated at the same time 
according to the state of their inputs and to their update 
function. Starting from some initial state, the network 
performs a trajectory in state space and eventually ar- 
rives on an attractor, where the same sequence of states 
is periodically repeated. Since the update rule is deter- 
ministic, the same state must always be followed by the 
same next state. If we represent the network states by 
points in the 2 7V -dimensional state space, each of these 
points has exactly one "output" , which is the successor 
state. We thus obtain a graph in state space. 

The size or length of an attractor is the number of 
different states on the attractor. The basin of attraction 
of an attractor is the set of all states that eventually 
end up on this attractor, including the attractor states 
themselves. The size of the basin of attraction is the 
number of states belonging to it. The graph of states 
in state space consists of unconnected components, each 
of them being a basin of attraction and containing an 
attractor, which is a loop in state space. The transient 
states are those that do not lie on an attractor. They are 
on trees leading to the attractors. 

Let us illustrate these concepts by studying the small 
K = 1 network shown in Figure 2.2, which consists of 4 
nodes: 




FIG. 2: A small network with K = 1 input per node. 



If we assign to the nodes 1,2,3,4 the functions invert, 
invert, copy, copy, an initial state 1111 evolves in the 
following way: 

1111 -> 0011 -> 0100 -> 1111 
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This is an attractor of period 3. If we interpret the bit se- 
quence characterizing the state of the network as a num- 
ber in binary notation, the sequence of states can also be 
written as 

15^3^4^ 15 
The entire state space is shown in Figure l2"Cl 

N S 
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FIG. 3: The state space of the network shown in Figure 2.2, if 
the functions copy, copy, invert, invert are assigned to the four 
nodes. The numbers in the squares represent states, and ar- 
rows indicate the successor of each state. States on attractors 
are shaded. 



There are 4 attractors, two of which 
(i.e., attractors of length 1). The sizes 
attraction of the 4 attractors are 6,6,2,2 
of node 1 is a constant function, fixing 
node at 1, the state of this node fixes 
network, and there is only one attractor, 
point. Its basin of attraction is of size 16. 
of the other nodes remain unchanged, 
then looks as shown in Figure 2.4. 

m H 

PC 9 ^ 

' r\- 



are fixed points 
of the basins of 
. If the function 
the value of the 
the rest of the 
which is a fixed 
If the functions 
the state space 




FIG. 4: The state space of the network shown in Figure 2.2, 
if the functions 1, copy, invert, invert are assigned to the four 
nodes. 



Before we continue, we have to make the definition of 
attractor more precise: as the name says, an attractor 
"attracts" states to itself. A periodic sequence of states 
(which we also call cycle) is an attractor if there are states 
outside the attractor that lead to it. However, some net- 
works contain cycles that cannot be reached from any 
state that is not part of it. For instance, if we removed 
node 4 from the network shown in Figure 2.2, the state 
space would only contain the cycles shown in Figure [2~Cl 
and not the 8 states leading to the cycles. In the follow- 
ing, we will use the word "cycle" whenever we cannot be 
confident that the cycle is an attractor. 



D. Applications 



Let us now make use of the definitions and concepts 
introduced in this section in order to derive some results 
concerning cycles in state space. First, we prove that in 
an ensemble of networks with update rule 1 (biased func- 
tions) or rule 2 (weighted classes), there is on an average 
exactly one fixed point per network. A fixed point is a 
cycle of length 1. The proof is slightly different for rule 
1 and rule 2. Let us first choose rule 2. We make use of 
the property that for every update function the inverted 
function has the same probability. The inverted function 
has all Is in the output replaced with Os, and vice versa. 
Let us choose a network state, and let us determine for 
which fraction of networks in the ensemble this state is a 
fixed point. We choose a network at random, prepare it 
in the chosen state, and perform one update step. The 
probability that node 1 remains in the same state after 
the update, is 1/2, because a network with the inverted 
function at node 1 occurs equally often. The same holds 
for all other nodes, so that the chosen state is a fixed 
point of a given network with probability 2~ N . This 
means that each of the 2 N states is a fixed point in the 
proportion 2~ N of all networks, and therefore the mean 
number of fixed points per network is 1. We will see 
later that fixed points may be highly clustered: a small 
proportion of all networks may have many fixed points, 
while the majority of networks have no fixed point. 

Next, we consider rule 1. We make now use of the 
property that for every update function a function with 
any permutation of the input states has the same proba- 
bility. This means that networks in which state A leads 
to state B after one update, and networks in which an- 
other state C leads to state B after one update, occur 
equally often in the ensemble. Let us choose a network 
state with n Is and -/V — n Os. The average number of 
states in a network leading to this state after one update 
is 2 N p n (l — p) N ~ n . Now, every state leads equally often 
to this state, and therefore this state is a fixed point in 
the proportion p n (l — p) N ~ n of all networks. Summation 
over all states gives the mean number of fixed points per 
network, which is 1. 

Finally, we derive a general expression for the mean 
number of cycles of length L in networks with K = 2 
inputs per node. The generalization to other values of K 
is straightforward. Let {Cl)n denote the mean number 
of cycles in state space of length L, averaged over the 
ensemble of networks of size N. On a cycle of length 
L, the state of each node goes through a sequence of 
Is and Os of period L. Let us number the 2 L possible 
sequences of period L of the state of a node by the index 
j, ranging from to to = 2 L — 1. Let rij denote the 
number of nodes that have the sequence j on a cycle of 
length L, and (Pl)\ & the probability that a node that 
has the input sequences I and k generates the output 
sequence j. This probability depends on the probability 
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distribution of update functions. Then 

!<■ ! ° i \ '> / 

(3) 

The factor 1/L occurs because any of the L states on the 
cycle could be the starting point. The sum is over all 
possibilities to choose the values {rij} such that J] ■ rij = 
N. The factor after the sum is the number of different 
ways in which the nodes can be divided into groups of the 
sizes no, 1%%, n 2 , . ■ . , n m . The product is the probability 
that each node with a sequence j is connected to nodes 
with the sequences I and k and has an update function 
that yields the output sequence j for the input sequences 
I and k. This formula was first given in the beautiful 
paper by Samuelsson and Troein [X 01 ] - 

We conclude this section with a picture of the state 
space of a network consisting of 10 nodes. 




FIG. 5: The state space of a network with 10 nodes 



E. Problems 

1. Show that the fraction 3/32 of all networks in the 
ensemble with N = 4 and K = 1 have the topology 
shown in Figure 2.2. 

2. Show that the fraction 3 3 /2 10 of all networks with 
the topology shown in Figure 2.2 have the state 
space topology shown in Figure 2.3, if the distri- 
bution of update functions is given by rule 1 with 
p=l/4. 

3. Which functions in Table [TT] correspond to the 
threshold functions in networks with K = 2, if we 
set h = ? 

4. Consider again K — 2 networks, and choose the up- 
date rules 3 (canalyzing functions), which are char- 
acterized by the parameters n, r, and p. Express 
the weight of each function in Tabic [TT] in terms of 
T], r, and p. 



5. Using Equation ([3]). show that in an ensemble of 
networks with update rule 1 or 2, there is on an 
average exactly one fixed point per network. 

3. ANNEALED APPROXIMATION AND 
PHASE DIAGRAMS 

The annealed approximation, which is due to Derrida 
and Pomeau [ll[ , is a useful tool to calculate certain net- 
work properties. It is a mean-field theory, which neglects 
possible correlations between nodes. The first assump- 
tion of the annealed approximation is that the network 
is infinitely large. This means that fluctuations of global 
quantities arc negligible. The second assumption of the 
annealed approximation is that the inputs of each node 
can be assigned at every time step anew. The following 
quantities can be evaluated by the annealed approxima- 
tion: 

1. The time evolution of the proportion of Is and 0s. 

2. The time evolution of the Hamming distance be- 
tween the states of two identical networks. 

3. The statistics of small perturbations. 

We will discuss these in the following in the order given in 
this list. One of the main results of these calculations will 
be the phase diagram, which indicates for which param- 
eter values the networks are frozen, critical or chaotic. 

A. The time evolution of the proportion of Is and 

Os 

Let bt denote the number of nodes in state 1, divided 
by TV. The proportion of nodes in state is then 1 — 
b t . We want to calculate 6 t+1 as function of b t within 
the annealed approximation. Since the K inputs of each 
node are newly assigned at each time step, the probability 
that to inputs of a node are in state 1 and the other 
inputs in state is — bt) K ~ m . Since we consider 

an infinitely large network, this probability is identical 
to the proportion of nodes that have m inputs in state 1. 

Let p m be the probability that the output value of a 
node with to inputs in state 1 is 1. Then we have 

b ^=J2 (^)pmbT(i-bt) K ~ m . (4) 

m=0 ^ ' 

If p m is independent of to, the right-hand side is iden- 
tical to p m , and b t reaches after one time step its sta- 
tionary value, which is the fixed point of Equation ([4"|. 
Among the above-listed update rules, this happens for 
rules 1 (biased functions) and 2 (weighted classes ) and 
4 (threshold functions). For rule 1, we have p m = 1/2, 
since the output values and 1 occur with equal proba- 
bility within each class of update functions. For rule 2, 
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we have p m = p by definition. For rule 4, the value of 
p m is independent of m because the value of Cjj is 1 and 
— 1 with equal probability, making each term Cij(2<jj — 1) 
to +1 and —1 with equal probability. Therefore p m is 
identical to the probability that the sum of K random 
numbers, each of which is +1 or —1 with probability 1/2, 
is at least as large as — h, 



K 



l>{K-h)/1 



Here, I is the number of +ls, and K — I the number of 
-Is. 

For rule 3 (canalyzing functions) we get Q 



h+i = h rjr + (1 - 6 t )(l - rj)r 
+b t (l - rj)p + (1 - b t )vp 
= r + r]{p — r) + b t (jp — r)(l 



2rj) . (5) 



The first two terms are the probability that the canalyz- 
ing input is on its canalyzing value, and that the output 
is then 1. The second two terms are the probability that 
the canalyzing input is not on its canalyzing value, and 
that the output is then 1. This is a one-dimensional map. 
The only fixed point of this map is 



b* = 



r + r}{p — r) 
l-(p-r)(l-277) 



Since the absolute value of the slope of this map is smaller 
than 1 everywhere, every iteration (|5|) will bring the value 
of bt closer to this fixed point. 

There exist also update rules where the fixed points 
are unstable and where periodic oscillations or chaos oc- 
cur. This occurs particularly easily when all nodes are 
assigned the same function (rule 5). For instance, if all 
nodes arc assigned the last one of the canalyzing func- 
tions occurring in the table of update functions [ill wc 
have the map 



h+i = l-b 2 t 



(0) 



The fixed point 



b* = 



-1 + V5 
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FIG. 6: The values of bt that still occur after the transient 
time for the map (3.5), as function of K. 



Let us consider if as a continous parameter. When it is 
increased, starting at 1, the map first has a stable fixed 
point and then shows a period-doubling cascade and the 
Feigenbaum route to chaos shown in Figure 13 Al [l2| . 

All these results for bt were derived within the annealed 
approximation, but they arc generally believed to apply 
also to the original networks with fixed connectivity pat- 
terns, if the thermodynamic limit is taken. If this is cor- 
rect, the following three statements are also correct: 

• All (apart from a vanishing proportion of) initial 
states with a given value of bo undergo the same 
trajectory b t with time. 

• This trajectory is the same for all networks (apart 
from a vanishing proportion) . 

• When time is so large that the dynamics have 
reached an attractor, the map bt+\{bt) is the same 
as in the initial stage for those values of b that can 
occur on the attractors. 

These assumptions appear plausible, since the paths 
through which a node can affect its own input nodes are 
infinitely long in a randomly wired, infinitely large net- 
work. Therefore we do not expect correlations between 
the update function assigned to a node and the states of 
its input nodes. Neither do we expect a correlation be- 
tween the function b t +i(b t ) and the question of whether 
a state is on an attractor. 



is unstable, since the slope of the map is (1 — y/E) at this 
fixed point, i.e., it has an absolute value larger than 1. 
The iteration ((6|) moves bt away from this fixed point, 
and eventually the network oscillates between all nodes 
being 1 and all nodes being 0. 

A map that allows for oscillations with larger period 
and for chaos is obtained for the update rule that the 
output is 1 only if all inputs are equal. This map is 
defined for general values of K and is given by 



UK 



(1 



b t ) K 



(7) 



B. The time evolution of the Hamming distance 

With the help of the Hamming distance, one can dis- 
tinguish between a frozen and a chaotic phase for RBNs. 
We make an identical copy of each network in the en- 
semble, and we prepare the two copies of a network in 
different initial states. The Hamming distance between 
the two networks is defined as the number of nodes that 
are in a different state. For the following, it is more con- 
venient to use the normalized Hamming distance, which 
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is the Hamming distance divided by N, i.e., the propor- 
tion of nodes that are in a different state, 

i=l 

If /it is very small, the probability that more than one 
input of a node differ in the two copies, can be neglected, 
and the change of h t during one time step is given by 

h t +i = Xht , (9) 

where A is called the sensitivity (l3| . It is K times the 
probability that the output of a node changes when one 
of its inputs changes. 

For the first four update rules listed in Section 2.2, the 
value of A is 

A = 2Kp(l — p) (biased functions) 

A = 1 — 8 (weighted classes, K = 1) 

A = a + 2/3 + 7= l + /3-<5 (weighted classes, K = 2) 

A = r(l-p) + (l-r)p+(K-l)(r](l-b t ) 

+ (1 — rf)bt)2p(l — p) (canalyzing functions) 

A = K J ^ (threshold functions) (10 

with / in the last line being the largest integer smaller 
than or equal to (K — h)/2. For rule 3 (canalyzing func- 
tions), the first two terms are the probability that the 
output changes when the canalyzing input is in a differ- 
ent state in the two network copies; the last term is the 
probability that the output changes when one of the other 
inputs is in a different state in the two copies, multiplied 
by the number of noncanalyzing inputs. This is the only 
one out of the 4 rules where the value of A depends on 6 t 
and therefore on time. 

The networks are in different phases for A < 1 and 
A > 1, with the critical line at A = 1 separating the two 
phases. In the following, we derive the properties of the 
networks in the two phases as far as possible within the 
annealed approximation. 

If A < 1, the normalized Hamming distance decreases 
to 0. If the states of the two copies differ initially in 
a small proportion of all nodes, they become identical 
for all nodes, apart from possibly a limited number of 
nodes, which together make a contribution to the nor- 
malized Hamming distance. A < 1 means also that if the 
two copies are initially in identical states and the state 
of one node in one copy is changed, this change propa- 
gates on an average to less than one other node. When 
the two copies differ initially in a larger proportion of 
their nodes, we can argue that their states also become 
identical after some time: we produce a large number Q 
of copies of the same network and prepare their initial 
states such that copy number q and copy number q + 1 
(for all q = 1,...,Q) differ only in a small proportion 
of their nodes. Then the states of copy number q and 



copy number q+l will become identical after some time, 
and therefore the states of all Q copies become identical 
(again apart from possibly a limited number of nodes). 
The final state at which all copies arrive must be a state 
where all nodes (apart from possibly a limited number) 
become frozen at a fixed value. If the final state was an 
attractor where a nonvanishing proportion of nodes go 
through a sequence of states with a period larger than 
1, different network copies could be in different phases 
of the attractor, and the normalized Hamming distance 
could not become zero. Ensembles with A < 1 are said 
to be in the frozen phase. 

All these considerations did not take into account that 
A itself may not be constant. For those update rules 
where bt assumes its fixed point value after the first time 
step, one can apply the reasoning of the previous para- 
graph starting at time step 2. For rule 3, the value bt 
approaches its fixed point more slowly, and therefore the 
value of A changes over a longer time period. It is there- 
fore possible that the Hamming distance shows initially 
another trend as during later times. Once bt has reached 
its fixed point value, A has become constant, and if then 
A < 1, the normalized Hamming distance will decrease 
to zero. In order to decide whether an ensemble is in 
the frozen phase, one must therefore evaluate A in the 
stationary state. 

For ensembles that have no stable stationary value of 
bt, the above considerations do not apply directly, since bt 
cycles through different values, and so does A. Further- 
more, the two copies may be in different phases of the 
cycle and will then never have a small normalized Ham- 
ming distance. For ensembles with a finite oscillation 
period T, one should evaluate the product of all values 
of A during one period. If this product is smaller than 1, 
a small normalized Hamming distance created in a copy 
of a network with a stationary oscillation, will decrease 
after one period. Using a similar reasoning as before, we 
conclude that then the normalized Hamming distance be- 
tween any two copies of the network will decrease to zero 
if they have initially the same value of b t . This means 
that all nodes (apart from possibly a limited number) go 
through a sequence of states that has the same period as 
b t . 

If A > 1 when bt has reached its stationary value, 
the normalized Hamming distance increases from then on 
with time and has a nonzero stationary value. A change 
in one node propagates on an average to more than one 
other node. If there is a fixed point or a short attractor, it 
is unstable under many possible perturbations. There is 
therefore no reason why all attractors should be short. In 
fact, attractors can be very long, and the ensemble is in 
a phase that is usually called chaotic, even though this 
is no real chaos because state space is finite and every 
trajectory becomes eventually periodic. When b t does 
not become stationary but periodic, we consider again 
the product of all values of A during one period. If this 
product is larger than 1, a small normalized Hamming 
distance between two copies with the same value of b t 
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will eventually become larger. This means that attrac- 
tors can be very long and need not have the period of 
h- 

For A = 1, the ensemble is at the boundary between 
the two phases: it is critical. A change in one node prop- 
agates on an average to one other node. The critical 
line can be obtained from Eqs. (fTU|) . leading to the phase 
diagram shown in Figure 13 Bl 
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FIG. 7: Phase diagram for the first 4 update rules (biased 
functions, weighted classes, canalyzing functions, threshold 
functions). Where there are more than 2 parameters, the 
remaining parameters were fixed at the values given in the 
respective graph titles. For theshold functions, the frozen 
phase is shaded. For K = 2, the model is critical between 
h = —2 and 2, for larger K, there exists no critical value, but 
the model is chaotic whenever it is not frozen. 



work. Let us consider again two identical networks, and 
let them be initially in the same state. Then let us flip 
the state of one node in the first network. One time step 
later, the nodes that receive input from this node differ in 
the two systems each with probability X/K = l/K (since 
A = 1 in a critical network). On an average, this is one 
node. Since the perturbation propagates to each node 
in the network with probability (K/N) * (X/K) = 1/N, 
the probability distribution is a Poisson distribution with 
mean value 1. We keep track of all nodes to which the 
perturbation propagates, until no new node becomes af- 
fected by it. We denote the total number of nodes af- 
fected by the perturbation by s. The size distribution of 
perturbations is a power law 



n(s) 



-3/2 



(13) 



for values of s that are so small that the finite system 
size is not yet felt, but large enough to see the power 
law. There are many ways to derive this power law. The 
annealed approximation consists in assuming that loops 
can be neglected, so that the perturbation propagates at 
every step through new bonds and to new nodes. In this 
case, there is no difference (from the point of view of the 
propagating perturbation) between a network where the 
connections are fixed and a network where the connec- 
tions are rewired at every time step. 

We begin our calculation with one "active" node at 
time 0, n a (t = 0) = 1, which is the node that is per- 
turbed. At each "time step" (which is different from real 
time!), we choose one active node and ask to how many 
nodes the perturbation propagates from this node in one 
step. These become active nodes, and the chosen node is 
now "inactive" . We therefore have a stochastic process 



All our results are based on a calculation for small ht- 
When h t is not infinitesimally small, ([9]) has the more 
general form 

ht+i = Xh t + uht + ■ • ■ , (11) 

with the highest power of ht being K, but we do not 
make here the effort to calculate the coefficent v or that 
of a higher-order term. We use this result only to obtain 
a relation between the stationary value of h t and the 
distance from the critical line: In the chaotic phase, but 
close to the critical line (where A is only slightly larger 
than 1), the stationary value of h t obtained from (fTTj) is 

h* = (X-l)/v. (12) 

It increases linearly with the distance from the critical 
line, as long as this distance is small. 

C. The statistics of small perturbations in critical 
networks 

Now let us have a closer look at the propagation of a 
perturbation that begins at one node in a critical net- 



n a (t+l)=n a (t) -1 + t 

for the number of "active" nodes, with £ being a random 
number with a Poisson distribution with mean value 1. 
The stochastic process is finished at time t = s when 
n a (s) = 0. s is the total number of nodes affected by the 
perturbation. 

Now we define Po(y,t) as the probability that the 
stochastic process has arrived at y = before or at time 
t, if it has started at n a — y at time t — 0. During the 
first step, y changes by Ay = £ — 1. If we denote the 
probability distribution of Ay with P(Ay), we obtain 

Po(y,t) = J d(Ay)P(Ay)P (y + Ay,t-l) 

~ J d(Ay)P(Ay) 

\ . & Pq I.. , 2 d 2 P dP ~ 
P o(y,t) + Ay^ + -(A y r^-— . 

The first term on the right-hand side cancels the left- 
hand side. The second term on the right-hand side is the 
mean value of Ay, which is zero, times d y Po. Were are 
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therefore left with the last two terms, which give after 
integration 



dP = i d 2 p 

dt 2 dy 2 



(14) 



This is a diffusion equation, and we have to apply the 
initial and boundary conditions 



Po(0,t) = 1 

Po(y,o) = o 

P (y, oo) = 1. 



(15) 



Expanding Pq in terms of eigenfunctions of the operator 
d/dt gives the general solution 

Po(y,t) =a + by+ / dwe'^ 4/4 (c w sin(wy) + d w cos(wy)) . 



The initial and boundary conditions fix the constants to 
a = 1 and gL = and c w = —2/ttlu. We therefore have 



1 



(16) 



which becomes for y = 1 



P (l,t) = l-~ du 

7T ./ W 



2 */4 



1 - G(t 



-1/21 



(17) 



for large i. The size distribution of perturbations is ob- 
tained by taking the derivative with respect to t, leading 
to (HH). 

Readers familar with percolation theory will notice 
that the spreading of a perturbation in a critical RBN 
is closely related to critical percolation on a Bcthe lat- 
tice. Only the probability distribution of the stochastic 
variable £ is different in this case. Since the result de- 
pends only on the existence of the second moment of y, 
it is not surprising that the size distribution of critical 
percolation clusters on the Bcthe lattice follows the same 
power law. 



D. Problems 

1. Explain why there is a finite critical region for K = 
2 and no critical value of A at all for K > 2 for 
update rule 4 (Figure [3B]) . 

2. For each of the 16 update functions for K = 2, con- 
sider an ensemble where all nodes are assigned this 
function. Find the function b t +i(bt). Find all fixed 
points of b and determine if they are stable. If bt be- 
comes constant for large times, determine whether 
the ensemble is frozen, critical or chaotic. If b t os- 
cillates for large times, determine whether the nor- 
malized Hamming distance between two identical 
networks that start with the same value of bo, goes 
to zero for large times. Interpret the result. 



3. If in a frozen network only a limited number of 
nodes may not be frozen for large times, and if in a 
chaotic network a nonvanishing proportion of nodes 
remain nonfrozen, what do you expect in a critical 
network? 



4. NETWORKS WITH K = 1 

Many properties of networks with K = 1 inputs per 
node can be derived analytically. Nevertheless, these 
networks are nontrivial and share many features with 
networks with larger values of K. Therefore it is very 
instructive to have a closer look at K = 1 networks. In 
this review, we will not reproduce mathematically exact 
results that require long calculations, as is for instance 
done in [3, [H| ■ Instead, we will present phenomenologi- 
cal arguments that reproduce correctly the main features 
of these networks and that help to understand how these 
features result from the network structure and update 
rules. We begin by studying the topology of K = 1 net- 
works. Then, we will investigate the dynamics on these 
networks in the frozen phase and at the critical point. Fi- 
nally, we will show that the topology of K = 1 networks 
can be mapped on the state space of K = N networks, 
which allows us to derive properties of the attractors of 
K = N networks, which are chaotic. 



A. Topology of K = 1 networks 

If each node has one input, the network consists of dif- 
ferent components, each of which has one loop and trees 
rooted in this loop, as shown in Figure 14 AL K = 1 net- 
works have the same structure as the state space pictures 
of other random Boolean networks, like the ones shown 
in Figures 2.3 and 2.4, only the arrows are inverted. 



FIG. 8: Example of a network with one input per node. It 
has two components, the larger component has a loop of size 
3 and two trees rooted in it (one of size 1 and one of size 6), 
and the smaller component has a loop of size 1 and one tree 
of size 1. 



Let us first calculate the size distribution of loops. We 
consider the ensemble of all networks of size N. In each 
network of the ensemble, each node chooses its input at 
random from all other nodes. The probability that a 
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given node is sitting on a loop of size I is therefore 



P(l) 



i 



N J N 



-2/N 



,-(l-l)/N 



N 



-l{l-l)/2N 

N 



-l 2 /2N 



(18) 



The first factor is the probability that the input to the 
first node is not this node. The second factor is the prob- 
ability that the input to the second node is not the first 
or second node, etc. The last factor is the probability 
that the input of the Ith node is the first node. The 
approximation in the second step becomes exact in the 
thermodynamic limit N — > oo for values of I that satisfy 
limjv— >oo l/N = 0. The approximation in the last step 
can be made if I is large. 

The probability that a given node is sitting on any loop 
is therefore proportional to 



P{l)dl ~ N 



-1/2 



This means that the total number of nodes sitting on 
loops is proportional to s/N. 

The mean number of nodes sitting on loops of size I in 
a network is 



NP(l) 



-l 2 /2N 



The cutoff in loop size is proportional to v N. For I <C 
y/~N, the mean number of nodes in loops of size I is 1. 
This result can also be obtained by a simple argument: 
The probability that a given node is sitting on a loop 
of size I is in the limit N — > oo simply 1/N, since the 
node almost certainly does not choose itself as input or 
as input of its input etc, but in the Ith step the first node 
must be chosen as input, which happens with probability 
1/N. 

The mean number of loops of size I in a network is 



NP(l) 
I 



,-l 2 /2N 



For I CvJVi this is simply l/l. Since loops are formed 
independently from each other in the limit TV — ► oo, the 
probability distribution of the number of loops of size I 
is a Poisson distribution with mean value l/l. 
The mean number of loops per network is 



J2np(i)/i 



-x 2 /2 



-dx 



1 



In TV 



means that the average tree size is proportional to vJV. 
The construction of a tree can be described formally ex- 
actly in the same way as we described the propagation 
of a perturbation in a critical network in the last section: 
we begin with a node sitting in a loop. The nodes that 
are not sitting in loops receive their input with equal 
probability from any node in the network. Our node 
is therefore chosen with probability 1/N by every node 
outside the loops as an input, and the probability distri- 
bution of the number of outputs into the tree is a Poisson 
distribution with mean value 1 (neglecting terms of the 
order TV" 1 / 2 ). In the same way, we find that the number 
of outputs of each of the newly found tree nodes is again 
a Poisson distribution with mean value 1. We iterate this 
process until we have identified all nodes that are part of 
this tree. 

The size distribution of trees is ~ s~ 3 / 2 . The cutoff 
must be s m ax ~ N in order to be consistent with what 
we know about the mean tree size and the total number 
of nodes in trees: The mean tree size is 



.S.N' 



" 3/2 ds ~ s 1/2 



N . 



The total number of nodes in trees is proportional to 



ss 



~ 3 ' 2 ds ~ N . 



B. Dynamics on K = 1 networks 

Knowing the topology of K = 1 networks, allows us to 
calculate their dynamical properties. After a transient 
time, the state of the nodes on the trees will be inde- 
pendent of their initial state. If a node on a tree does 
not have a constant function, its state is determined by 
the state of its input node at the previous time step. All 
nodes that are downstream of a node with a constant 
function will become frozen. If there is no constant func- 
tion in the loop and the path from the loop to a node, 
the dynamics of this node is slaved to the dynamics of 
the loop. 

If the weight of constant functions, 6, is nonzero, the 
probability that a loop of size I does not contain a frozen 
function is (1 — 6) , which goes to zero when I is much 
larger than 1/6. Therefore only loops smaller than a 
cutoff size can have nontrivial dynamics. 

The number and length of the attractors of the net- 
work are determined by the nonfrozen loops only. Once 
the cycles that exist on each of the nonfrozen loops are 
determined, the attractors of the entire networks can be 
found from combinatorial arguments. 



for large N. This is identical to the mean number of 
components. 

Next, let us consider the trees rooted in the loops. 
There are of the order of N nodes, which sit in oc \/N 
trees, each of which is rooted in a relevant node. This 



1. Cycles on loops 

Let us therefore focus on a loop that has no constant 
function. If the number of "invert" functions is odd, we 
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call the loop an odd loop. Otherwise it is an even loop. 
Replacing two "invert" functions with copy functions and 
replacing the states (Ti(t) of the two nodes controlled by 
these functions and of all nodes in between with 1 — Oi(£), 
is a bijective mapping from one loop to another. In par- 
ticular, the number and length of cycles on the loop is 
not changed. All odd loops can thus be mapped on loops 
with only one "invert" function, and all even loops can 
be mapped on loops with only "copy" functions. 

We first consider even loops with only "copy" func- 
tions. These loops have two fixed points, where all nodes 
are in the same state. If I is a prime number, all other 
states belong to cycles of period I. Any initial state oc- 
curs again after I time steps. Therefore the number of 
cycles on an even loop is 



(19) 



/ 



if I is a prime number. The numerator counts the num- 
ber of states that are not fixed points. The first term is 
therefore the number of cycles of length /. Adding the 
two fixed points gives the total number of cycles. If I is 
not a prime number, there exist cycles with all periods 
that are a divisor of I. 

Next, let us consider odd loops with one "invert" func- 
tion. After 21 time steps, the loop is in its original state. 
If I is a prime number, there is only one cycle that has a 
shorter period. It is a cycle with period 2, where at each 
site Os and Is alternate. The total number of cycles on 
an odd loop with a prime number / is therefore 



21 



1 



(20) 



If I is not a prime number, there are also cycles with a 
period that is twice a divisor of I. 



2. K = 1 networks in the frozen phase 



For networks with K 
eter A is 



1 input per node, the param- 



A = 1 - 8. 



(21) 



Thefore, only networks without constant functions are 
critical. Networks with 8 > are in the frozen phase. 
The mean number of nonfrozen nodes on nonfrozen loops 
is given by the sum 



1-8 



(22) 



We call these loops the relevant loops. We call the nodes 
on the relevant loops the relevant nodes, and we denote 
their number with N re i. The mean number of relevant 
loops is given by the sum 



1-*)' ~ln£- 



(23) 



with the last step being valid for small 8. 

The probability that the activity moves up the tree to 
the next node is 1 — 8 at each step. The mean number of 
nonfrozen nodes on trees is therefore 



(24) 



and the total mean number of nonfrozen nodes is (1 — 
8)/8 2 . This is a finite number, which diverges as 5~ 2 
when the critical point 8 = is approached. 



3. Critical K = 1 networks 

If the proportion of constant functions 8 is zero, the 
network is critical, and all loops arc relevant loops. There 
are no nodes that are frozen on the same value on all 
attractors. A loop of size 1 has a state that is constant 
in time, but in can take two different values. Larger loops 
have also two fixed points, if they are even. Part of the 
nodes in a critical K = 1 networks are therefore frozen on 
some attractors or even on all attractors, however, they 
can be frozen in different states. 

The network consists of ~ ln7V/2 loops, each of which 
has of the order 2 l /I cycles of a length of the order I. The 
size of the largest loop is of the order of y/N. The num- 
ber of attractors of the network results from the number 
of cycles on the loops. It is at least as large as the prod- 
uct of all the cycle numbers of all the loops. If a cycle 
is not a fixed point, there are several options to choose 
its phase, and the number of attractors of the network 
becomes larger than the product of the cycle numbers. 
An upper bound is the total number of states of all the 



loops, which is 2 



N r . 



and a lower bound is the 



number of attractors on the largest loop, which is of the 
order e bVW /VN > e b '^ with b' < b < a. From this 
it follows that the mean number of attractors of critical 
K = 1 networks increases exponentially with the number 
of relevant nodes. A complementary result for the num- 
ber of cycles (Cr) of length L, which is valid for fixed L 
in the limit N — ► oo is obtained by the following quick 
calculation: 



(C 



L/N 




Kl c 



,{H L -l)\nVN _ N (H L -l)/2 



(25) 



Here, rij is the number of loops of size I, l c is the cutoff 
in loop size oc \/~N, and ki is the number of states on a 
loop of size I that belong to a cycle of length L. This 
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is zero for many loops. The average over an /-interval of 
size L is identical to Hi,, which is the number of cycles 
on an even loop of size L. A more precise derivation of 
this relation, starting from the K = 1 version of ^ and 
evaluating it by making a saddle-point approximation, 
can be found in [l||, which is inspired by the equivalent 
calculation for K = 2 critical networks in [Tc| . 

The length of an attractor of the network is the least 
common multiple of the cycle lengths of all the loops. A 
quick estimate gives 

jya log N 

since the length of the larger loops is proportional to 
V~N, and this has to be taken to a power which is the 
number of loops. A more precise calculation [l7T ] gives 
this expression, multiplied with a factor N b / log N, which 
does not modify the leading dependence on N. 



for the dependence on N of the number and size of the 
basins of attraction of the different attractors. 

Let us first consider networks in the frozen phase. As 
we have seen, there is at most a limited number of small 
nonfrozen loops. Their number is independent of system 
size, and therefore the number of attractors is also inde- 
pendent of the system size. The initial state of the nodes 
on these nonfrozen loops determines the attractor. The 
initial states of all other nodes are completely irrelevant 
at determining the attractor. 

The size of the basin of attraction of an attractor is 
therefore 2 N ~ Nr " l 7 multiplied with the length of the at- 
tractor, i.e., it is 2 N , divided by a factor that is indepen- 
dent of N. The proportion of state space belonging to a 
basin is therefore also independent of N. If we define the 
basin entropy [HI by 



S = -y^Pq \np a 



(26) 



C. Dynamics on K — N networks 

The topology of a K = 1 network is identical to the 
topology of the state space of a K = N network, when 
all update functions are chosen with the same weight. 
The reason is that in a K = N network, the state that 
succeeds a given state can be every state with the same 
probability. Thus, each state has one successor, which is 
chosen at random among all states. In the same way, in 
a K = 1 network, each node has one input node, which 
is chosen at random among all nodes. The state space 
of a K = N network consists of 2 N nodes, each of which 
has one successor. We can now take over all results for 
the topology of K = 1 networks and translate them into 
state space: 

The K = N networks have of the order of log(2 7V ) oc N 
attractors. The largest attractor has a length of the order 
V2 N = 2 N / 2 , and this is proportional to the total number 
of states on attractors. All other states are transient 
states. An attractor of length I occurs with probability 

1/1 in<2 w / 2 . 

Clearly, K = N networks, where all update functions 
are chosen with the same weight, arc in the chaotic phase. 
The mean number of nodes to which a perturbation of 
one node propagates, is N/2. At each time step, half the 
nodes change their state, implying also that the network 
is not frozen. 



D. Application: Basins of attraction in frozen, 
critical and chaotic networks 

The advantage of K = 1 networks is that they are an- 
alytically tractable and can teach us at the same time 
about frozen, chaotic and critical behavior. We will dis- 
cuss in the next section to what extent the results apply 
to networks with other values of K. Based on our in- 
sights into K = 1 networks, we derive now expressions 



with p a being the fraction of state space occupied by the 
basin of attraction of attractor a, we obtain 

S = const 

for a K = 1 network in the frozen phase. 

Next, let us consider the chaotic K = N network en- 
semble. There are on an average l/l attractors of length 
I, with a cutoff around 2 N I 2 . The basin size of an attrac- 
tor of length / is of the order 12 N / 2 , which is / times the 
average tree size. The basin entropy is therefore 



S 



1 



Z-^i i 2 N / 2 2 N I 2 



log; 



/ 



2-N/2 



log xdx = const 



(27) 

Finally, we evaluate the basin entropy for a critical 



K = 1 network. There arc of the order e 



aVN 



attractors 



with approximately equal basin sizes, and therefore the 
basin entropy is 



S 



N CX N re l 



(28) 



While frozen and chaotic networks have a finite basin 
entropy, the basin entropy of critical networks increases 
as the number of relevant nodes LL8l . 



E. Problems 

1. How many cycles does an even (odd) loop of size 6 
have? 

2. Count the attractors of the network shown in Fig- 
ure 4.1 for all four cases where loop 1 and/or loop 
2 are even/odd. 

3. How does the transient time (i.e. the number of 
time steps until the network reaches an attrac- 
tor) increase with N for (a) K = 1 networks in 
the frozen phase, (b) critical K — 1 networks, (c) 
chaotic K — N networks? 
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4. Consider the subcnscmble of all critical K = 1 
networks that have the same wiring, but all pos- 
sible assignments of "copy" and "invert" functions. 
Which property determines the probability that a 
network has a fixed point attractor? If it has such 
an attractor, how many fixed point attractors does 
the network have in total? Conclude that there is 
on an average one fixed point per network in this 
subensemble. 

5. Verify the identity k\ = Hl used in calculation (|2"5l) . 

6. How does the basin entropy for K = 1 networks 
depend on the parameter 5 when 5 becomes very 
small? Find an answer without performing any cal- 
culations. 



5. CRITICAL NETWORKS WITH K = 2 

In the previous section, we have derived many proper- 
ties of frozen, critical and chaotic networks by studying 
ensembles with K = 1. Many results are also valid for 
RBNs with general values of K. In this section, we focus 
on critical K = 2 networks. These networks, as well as 
critical networks with larger values of K, differ in one im- 
portant respect from critical K = 1 networks: they have 
a frozen core, consisting of nodes that are frozen on the 
same value on all attractors. We have obtained this result 
already with the annealed approximation: The normal- 
ized Hamming distance between two identical networks 
is close to the critical point given by (fl"2|) . which means 
that it is zero exactly at the critical point. For K = 1, 
there exists no chaotic phase and no Equation (Tj"2"j) . and 
therefore the observation that all nodes may be nonfrozen 
in critical K = 1 networks is not in contradiction with 
the annealed approximation. 

We will first explain phcnomcnologically the features 
of critical K = 2 networks, and then we will derive some 
of these features analytically. 

A. Frozen and relevant nodes 

The frozen core arises because there are constant func- 
tions that fix the values of some nodes, which in turn 
lead to the fixation of the values of some other nodes, 
etc. Let us consider Figure 15 Al as an example. This 
network has the same number of constant and reversible 
functions, as is required for critical networks (although 
this classification only makes sense for large networks, 
where the thermodynamic limit becomes visible). Node 
5 has a constant function and is therefore frozen on the 
value (indicated by a darker grey shade) after the first 
time step. Node 6 has a canalyzing function which gives 
1 as soon as one of the inputs is 0. Therefore node 6 is 
frozen in state 1 (indicated by a lighter grey shade) no 
later than after the second time step. Then node 7 has 
two frozen inputs and becomes therefore also frozen. Its 
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FIG. 9: Example of a network with 8 nodes and K = 2. 
The functions are those of Table 2.2 (but written horizontally 
instead of vertically), with the first input node being the one 
with the lower number. 



value is then 1. Node 8 has a canalyzing function which 
gives as soon as one of the inputs is 1, and will there- 
fore end up in state 0. These four nodes constitute the 
frozen core of this network. At most after 4 time steps, 
each of these nodes assumes its stationary value. If we 
remove the frozen core, we are left with a K = 1 network 
consisting of nodes 1 to 4, with "copy" and "invert" func- 
tions between these nodes. For instance, node 4 copies 
the state of node 2 if node 7 is in state 1. Node 3 copies 
the state of node 2, node 1 inverts the input it receives 
from node 3, and node 2 inverts the input it receives from 
node 1. The nodes 1,2,3 form an even loop, and node 4 is 
slaved to this loop. Nodes 1,2,3 are therefore the relevant 
nodes that determine the attractors. We can conclude 
that this network has 4 attractors: two fixed points and 
two cycles of length 3. 

There is a different mechanism by which a frozen core 
can arise, which is illustrated by assigning another set 
of update functions to the same network, as shown in 
Figure 5.2. This network contains only canalizing update 

1 101 
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FIG. 10: A network with the same topology as the previous 
network, but with only canalyzing functions. 

functions of the type C2, and such a network could be 
classified as critical if it was much larger. We begin again 
by fixing node 5 at value 1, and we denote this as 5i. In 
the next time step, this node may have changed its state, 
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but then node 7 will be in state 1, because it is canalyzcd 
to this value by node 5. By continuing this consideration, 
we arrive at the following chain of states: 

5i -» 7i -» 4i -> 5 . 

This means that node 5 must eventually assume the state 
0, and we continue from here by following again canalyz- 
ing connections: 

5o - 6i -> 7i -> (4i,8o) -> (5 0j 6i) -> (6i,7i) 
(4 l5 7i, 8 ) -► (4i, 5 , 6i, 8 ) -> (5 , 61,70 

-> (4i,6i,7i,8 ) (4 1 ,5 ,6i,7 1 ,8o) (4 X , 5 , 6i, 7 U 8 ) 

From this moment on, nodes 4 to 8 are frozen. Nodes 
1,2,3 form a relevant loop with the functions invert, in- 
vert, copy, just as in the previous example. 

In order to better understand how the frozen core arises 
in this case, consider the loop formed by the nodes 6,7,8: 
This is a self-freezing loop. If the nodes 6,7,8 are in the 
states 1,1,0, they remain forever in these states, because 
each node is canalyzcd to this value by the input it re- 
ceives within the loop. This loop has the same effect on 
the network as have nodes with constant functions. Once 
this loop is frozen, nodes 4 and 5 become also frozen. One 
can imagine networks where such a loop never freezes, 
but this becomes very unlikely for large networks. 

The networks shown in the previous two figures were 
designed to display the desired properties. In general, 
small networks differ a lot in the number of frozen and 
nonfrozen nodes, as well as in the size and structure of 
their relevant component(s) and their attractors. The 
specific properties particular to the frozen and chaotic 
phase and to the critical line become clearly visible only 
for very large networks. 

A network of intermediate size is the basis of Figure 
5.3, which shows the nonfrozen part of a critical K = 2 
network with 1000 nodes. There are 100 nonfrozen nodes 
in this network, indicating that the majority of nodes are 
frozen. Among the 100 nonfrozen nodes, only 5 nodes are 
relevant, and only 6 nodes have two nonfrozen inputs. 
The relevant nodes are arranged in 2 relevant compo- 
nents. They determine the attractors of the network, 
while all other nodes sit on outgoing trees and are slaved 
to the dynamics of the relevant nodes. This figure resem- 
bles a lot a K = 1 network. The only difference is that 
there are a few nodes with two inputs. 

Analytical calculations, part of which are explained in 
the next section, give the following general results for 
critical K = 2 networks in the thermodynamic limit N — > 
oo: 

1. The number of nodes that do not belong to the 
frozen core, is proportional to TV 2 / 3 for large N. 

2. If the proportion of nodes with a constant function 
is nonzero, the frozen core can be determined by 
starting from the nodes with constant functions and 
following the cascade of freezing events. 




FIG. 11: The nonfrozen part of a K = 2 network with 1000 
nodes. Shown are the 100 nonfrozen nodes. 5 nodes are rele- 
vant (white), and 6 nodes (black) have two nonfrozen inputs. 



3. If the proportion of nodes with a constant function 
is zero (which means that the network contains only 
canalyzing functions), the frozen core can be deter- 
mined by starting from self-freezing loops. 

4. The number of nodes that are nonfrozen and that 
receive 2 nonfrozen inputs is proportional to iV 1 / 3 . 

5. The number of relevant nodes is proportional to 
iV 1 / 3 . They are connected to relevant components, 
which consist of loops and possibly additional links 
within and between the loops. 

6. The number of relevant nodes that have two rele- 
vant inputs remains finite in the limit N — ► oo. 

7. The number of relevant components increases as 
log TV 1 / 3 . 

8. The cutoff of the size of relevant components scales 
as N 1 / 3 . 

The complete list of these results is given in 1 1911 . but part 
of the results can be found in earlier papers [20, [2l|, [23] ■ 



B. Analytical calculations 

After this qualitative introduction to critical networks, 
let us derive the main results for the scaling of the num- 
ber of nonfrozen and relevant nodes with N. Computer 
simulations of critical networks show the true asymptotic 
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scaling only for very larger networks with more than 
100000 nodes. For this reason, the values 2/3 and 1/3 
for the critical exponents characterizing the number of 
nonfrozen and relevant nodes has been known only since 
2003. 

Flyvbjerg [23[ was the first one to use a dynamical 
process that starts from the nodes with constant update 
functions and determines iteratively the frozen core. Per- 
forming a mean- field calculation for this process, he could 
identify the critical point. We will go now beyond mean- 
field theory. 

We consider the ensemble of all K = 2 networks of 
size N with update rule 2 (weighted functions), where 
the weights of the Ci, reversible, C2 and constant func- 
tions are a, /3, 7 and S. These networks are critical for 
[3 = 5. We begin by assigning update functions to all 
nodes and by placing these nodes according to their func- 
tions in four containers labelled J-, C\, C2, and TZ. These 
containers then contain Nf, N Cl , N C2 , and N r nodes. We 
treat the nodes in container C\ as nodes with only one 
input and with the update functions "copy" or "invert" . 
As we determine the frozen core, the contents of the con- 
tainers will change with time. The "time" we are defining 
here is not the real time for the dynamics of the system. 
Instead, it is the time scale for the process that we use 
to determine the frozen core. One "time step" consists 
in choosing one node from the container T, in selecting 
the nodes to which this node is an input, and in deter- 
mining its effect on these nodes. These nodes change 
containers accordingly. Then the frozen node need not 
be considered any more and is removed from the system. 
The containers now contain together one node less than 
before. This means that container T contains only those 
frozen nodes, the effect of which on the network has not 
yet been evaluated. The other containers contain those 
nodes that have not (yet) been identified as frozen. The 
process ends when container T is empty (in which case 
the remaining nodes are the nonfrozen nodes), or when 
all the other containers are empty (in which case the en- 
tire network freezes). The latter case means that the 
dynamics of the network go to the same fixed point for 
all initial conditions. 

This process is put into the following equations, which 
describe the changes of the container contents during one 
"time step". 



AN r = 

&N C2 = 
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an input with probability 2/N and becomes then a C\- 
node. This explains the first equation and the first term 
in the third equation. Each node in container C2 chooses 
the selected frozen node as an input with probability 
2/N. With probability 1/2, it then becomes frozen, be- 
cause the frozen node is with probability 1 /2 in the state 
that fixes the output of a C2-node. If the C2-node does 
not become frozen, it becomes a Ci-node. This explains 
the terms proportional to N C2 . Each node in container C\ 
chooses the selected frozen node as an input with prob- 
ability 1/N. It then becomes a frozen node. Finally, 
the —1 in the equation for ANf means that the chosen 
frozen node is removed from the system. In summary, the 
total number of nodes, N, decreases by one during one 
time step, since we remove one node from container T . 
The random variable £ captures the fluctuations around 
the mean change ANf. It has zero mean and variance 
(JV C1 +iV C2 )/iV. The first three equations should contain 
similar noise terms, but since the final number of nodes of 
each class is large for large N, the noise can be neglected 
in these equations. We shall see below that at the end 
of the process most of the remaining nodes are in con- 
tainer C± , with the proportion of nodes left in containers 
C2 and TZ vanishing in the thermodynamic limit. Figure 
] illustrates the process of determining the frozen core. 
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The terms in these equations mean the following: Each 
node in container TZ chooses the selected frozen node as 



FIG. 12: Illustration of the freezing process. (1) Initially, a 
frozen node is chosen (marked in white), (2) then it is deter- 
mined to which node(s) this is an input and the effect on those 
nodes is determined. (3) Then, the selected frozen node is re- 
moved. (4) The last picture sketches the final state, where all 
frozen nodes have been removed and most remaining nodes 
have 1 nonfrozen input. 



The number of nodes in the containers, N, can be used 
instead of the time variable, since it decreases by one 
during each step. The equations for N r and N C2 can 
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then be solved by going from a difference equation to a 
differential equation, 



AN r ^ 
AN ~ 

which has the solution 

(3N 2 



dN r 



2N r 



N r 



N' 



1 N 2 



(30) 



where we have now denoted the total number of nodes 
with N mi , since the value of N changes during the pro- 
cess. Similarly, we find if we neglect the noise term for a 
moment 



Nf = N(8-j3) + 



(3N 2 



N C1 = N(a + 7 + 2(3) - 2 



N 2 {P + l) 



(31) 



From this result, one can derive again the phase dia- 
gram, as we did by using the annealed approximation: 
For 5 < (3, i.e. if there are more frozen than reversible 
update functions in the network, we obtain Nf = at a 
nonzero value of N, and the number of nonfrozen nodes 
is proportional to N ml . We are in the chaotic phase. For 
S > (3, there exists no solution with Nf = and N > 0. 
The network is in the frozen phase. For the critical net- 
works that we want to focus on, we have 8 = (3, and the 
process stops at Nf = 1 = (3N 2 /N mi if we neglect noise. 
This means that N — ^jN in% / j3 at the end of the pro- 
cess. The number of nonfrozen nodes would scale with 
the square root of the network size. This is not what is 
found in numerical studies of sufficiently large networks. 
We therefore must include the noise term. Noise becomes 
important only after Nf has become small, when most 
nodes are found in container C\, and when the variance 
of the noise has become unity, (£ 2 ) = 1. Inserting the 
solution for N r into the equation for Nf, we obtain then 



dN 
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with the step size dN = 1. We want to transform this into 
a Fokkcr-Planck-cquation. Let P(Nf,N) be the prob- 
ability that there are Nf nodes in container T at the 
moment where there are N nodes in total in the contain- 
ers. This probability depends on the initial node number 
Ni n i, and on the parameter (3. The sum 

oo poo 

]T P(N f ,N)~ / P(Nf,N)dN f 



N 



is the probability that the stochastic process is not yet 
finished, i.e. the probability that Nf has not yet reached 
the value at the moment where the total number of 
nodes in the containers has decreased to the value N. 
This means that systems that have reached Nf = must 
be removed from the ensemble, and we therefore have to 



impose the absorbing boundary condition P(Q,N) = 0. 
Exactly in the same way as with calculation (114jl . we 
obtain then 

dP _ 8 (Nf , 13N\ ld 2 P 

F+ 28N 2 ' {M) 



dN dN f \ N N 
We introduce the variables 



Nf N 
x = —?= an d V = /„ T - - , ^o/? (34) 
v/jV (N mt /0) 2 / 3 

and the function f(x, y) = (N ini /j3)^ s P(N f ,N). We 
will see in a moment that f(x,y) does not depend ex- 
plicitely on the parameters N ml and (3 with this defini- 
tion. The Fokker-Planck equation then becomes 

*f+/+(t+* a/a )S+^=- <»> 

Let W(N) denote the probability that N nodes are left 
at the moment where Nf reaches the value zero. It is 

/•oo /• oo 

W{N) = I P(Nf,N)dN f - I P(Nf,N-l)dN f . 
Jo Jo 

Consequently, 

W(N) = A J°° P{Nf,N)dN f 

8 f°° 
= (N ini /f3)-V 3 —VN J f(x,y)dx 

= (N ini /l3)- 2 ^-^y f(x,y)dx 



{N ini /(3)- 2 ' s G{y) 



(36) 



with a scaling function G(y). W(N) must be a normal- 
ized function, 



W(N)dN 



G(y)dy = 1 



This condition is independent of the parameters of the 
model, and therefore G(y) and f(x, y) are independent 
of them, too, which justifies our choice of the prefactor in 
the definition of /(a;, y). The mean number of nonfrozen 
nodes is therefore 



NW(N)dN = (N ini //3) 2/3 / G{y)ydy 



N 

Jo JO 

(37) 

which is proportional to (N ml / f3) 2 ' 3 . From Equations 
POP and the corresponding equation for the C2-nodes we 
find then that the number of nonfrozen nodes with two 
nonfrozen inputs is proportional to N 1 ^ 3 . This is a van- 
ishing proportion of all nonfrozen nodes. 

The nonfrozen nodes receive their (remaining) input 
from each other, and we obtain the nonfrozen part of 
the network by randomly making the remaining connec- 
tions. If we neglect for a moment the second input of 
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those nonfrozen nodes that have two nonfrozcn inputs, 
we obtain a K = 1 network. The number of relevant 
nodes must therefore be proportional to the square root 
of number of nonfrozen nodes, i.e. it is N re i ~ TV 1 / 3 , 
and the number of relevant components is of the order 
In A 1 / 3 , with the largest component of the order of A 2 / 3 
nodes (including the trees). Adding the second input to 
the nonfrozen nodes with two nonfrozen inputs docs not 
change much: The total number of relevant nodes that 
receive a second input is a constant (since each of ~ TV 1 / 3 
relevant nodes receives a second input with a probabil- 
ity proportional to A -1 / 3 ). Only the largest loops are 
likely to be affected, and therefore only the large relevant 
components may have a structure that is more complex 
than a simple loop. Most nonfrozen nodes with two non- 
frozen inputs sit in the trees, as we have seen in Figure 
5.3. The mean number and length of attractors can now 
be estimated in the following way: The attractor num- 
ber must be at least as large as the number of cycles on 
the largest relevant loop, and therefore it increases expo- 
nentially with the number of relevant nodes. The mean 
attractor length becomes larger as for K = 1 networks, 
since complex relevant components can have attractors 
that comprise a large part of their state space, as was 
shown in [24|. Such components arise with a nonvanish- 
ing probability, and they dominate therefore the mean 
attractor length, which therefore increases now exponen- 
tially with the number of relevant nodes. 

The conclusions derived in the last paragraph can be 
made more precise. Interested readers arc referred to 

All these results are also valid for K = 2 networks 
with only canalyzing functions. As mentioned before, 
the frozen core of canalyzing networks arises through self- 
freezing loops. The resulting power laws are the same as 
for networks with constant functions, as was shown in 

C. Problems 

1. What is the number of attractors of the network 
shown in Figure 5.3 for all four cases where loop 1 
and/or loop 2 are even/odd? 

2. Assume there are 4 relevant nodes, one of them 
with two relevant inputs. List all topologically dif- 
ferent possibilities for the relevant components. 

3. Using Equation ([36]) . figure out how the probability 
that the entire network freezes depends on A. 

6. NETWORKS WITH LARGER K 

Just as we did for K = 2, we consider larger values of 
K only for those update rules that lead to fixed points of 
bt (i.e. of the proportion of Is), and therefore to a critical 
line separating a frozen and a chaotic phase. 



Let us first consider the frozen phase, where the sensi- 
tivity A is smaller than 1. The probability that a certain 
node is part of a relevant loop of size I is for large A 
obtained by the following calculation: the node has K 
inputs, which have again each K inputs, etc., so that 
there are K 1 ^ 1 nodes that might choose the first node as 
one of its K inputs, leading to a connection loop. The 
chosen node is therefore part of K l /N connection loops 
of length I on an average. The probability that a given 
connection loop has no frozen connection is (X/K) 1 , and 
therefore the mean number of relevant loops of size I is 
A /I. The mean number of relevant nodes is then 

<AW=E Ai = 7^y ( 38 ) 
i 

This is the same result as (|2"2"|) . which we derived for 
K = 1. The mean number of nonrelevant nodes to which 
a change of the state of a relevant node propagates is 
given by the same sum, since in each step the change 
propagates on an average to A nodes. By adding the 
numbers of relevant and nonrelevant nonfrozen nodes, 
we therefore obtain again a mean number of A/(l — A) 2 
nonfrozen nodes, just as in the case K = 1. We conclude 
that the frozen phases of all RBNs are very similar. 

Now we consider critical networks with K > 2. The 
number of nonfrozen nodes scales again as N 2 ' 3 and the 
number of relevant nodes as A 1 / 3 . The number of non- 
frozen nodes with k nonfrozen inputs scales with A as 
N (3-k)/3_ Thege resultg 

are obtained by generalizing the 
procedure used in the previous section for determining 
the frozen core [US- By repeating the considerations of 
the previous paragraph with the value A = 1, we find 
that in all critical networks the mean number of relevant 
loops of size i is l/l - as long as / is smaller than a cutoff, 
the value of which depends on A. For K = 1 the cutoff is 
at y/N, for K = 2, it is at A 1 / 3 , and this value does not 
change for larger K. There exists a nice phcnomcnologi- 
cal argument to derive the scaling ~ N 2 / 3 of the number 
of nonfrozen nodes [27| : The number of nonfrozen nodes 
should scale in the same way as the size of the largest per- 
turbation, since the largest perturbation affects all nodes 
on the largest nonfrozen component. The cutoff s max in 
the size of perturbations (see Equation (fT3|) ) is given by 
the condition that n(s max ) ~ 1/N. Perturbations larger 
than this size occur only rarely in networks of size A, 
since n(s) is the probability that a perturbation of one 
specific node (out of A the nodes) affects s nodes in total. 
Using Equation (fl"3f . we therefore obtain 

Smax ~ A 2 / 3 . (39) 

This argument does not work for K = 1, where we have 
obtained s max ~ A in section 01 The reason is that 
critical networks with K = 1 have no frozen core, but 
every node that receives its input from a perturbed node 
will also be perturbed. 

As far as the chaotic phase is concerned, there are good 
reasons to assume that it displays similar features for all 
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K. We have explicitely considered the case K = TV . Nu- 
merical studies show that the basin entropy approaches 
a constant with increasing K also when the value of K is 
fixed [l|[ . When A is close to 1 , there is a frozen core that 
comprises a considerable part of the network. We can ex- 
pect that the nonfrozen part has a state space structure 
similar to that of the K = TV networks. 



interesting question to address in this context is whether 
the networks remain in the neighborhood of one attrac- 
tor (which can be tested by evaluating the return prob- 
ability after switching off the noise), or whether they 
move through large regions of state space. Investigations 
of networks with such a type of noise can be found in 




7. OUTLOOK 

There are many possibilities of how to go beyond RBNs 
with synchronous update. In this last section, we will 
briefly discuss some of these directions. 

A. Noise 

Synchronous update is unrealistic since networks do 
not usually have a central pacemaker that tells all nodes 
when to perform the next update. Asynchronous up- 
date can be done either deterministically by assigning to 
each node an update time interval and an initial phase 
(i.e. the time until the first update), or stochastically 
by assigning to each node a time-dependent probability 
for being updated. We focus here on stochastic update, 
since all physical systems contain some degree of noise. 
In particular, noise is ubiquitous in gene regulatory net- 
works [28| . Boolean networks with stochastic update are 
for instance investigated in [2^, [3(| • The frozen core ob- 
viously remains a frozen core under stochastic update, 
and the relevant nodes remain relevant nodes. The most 
fundamental change that occurs when one switches from 
deterministic to stochastic update is that there is now in 
general more than one successor to a state. The set of 
recurrent states comprises those states that can reoccur 
infinitely often after they have occurred for the first time. 
However, if there is a path in state space from each state 
to a fixed point or to a cycle that has only one succes- 
sor for each state, the network behaves deterministically 
for large times, in spite of the stochastic update. This 
occurs in networks where all relevant nodes sit on loops: 
an even loop has two fixed points, and an odd loop has 
an attractor cycle of length 21, where each state has only 
one successor in state space (apart from itself). If the 
number of relevant loops increases logarithmically with 
system size TV, the number of attractors then increases 
as a power law of TV. This means that critical K = 1 net- 
works with asynchronous update have attractor numbers 
that increases as a power law with system size. In (30| 
it is argued that in critical K = 2 networks, where not 
all relevant components are simple loops, the attractor 
number is still a power law in TV. 

The situation becomes different when the noise does 
not only affect the update time but also the update func- 
tion. Then the output of a node can deviate from the 
value prescribed by the update function with a proba- 
bility that depends on the strength of the noise. The 



B. Scale-free networks and other realistic network 
structures 

Real networks do not have a fixed number of inputs per 
node, but do often have a power-law distribution in the 
number of inputs or the number of outputs [33j | . Boolean 
dynamics on such networks has been studied 34], how- 
ever, how this affects the power laws in critical networks, 
is only partially known [27] . 

There are many more characteristics of real networks 
that are not found in random network topologies, such 
as clustering, modularity, or scale invariance. The effect 
of all these features on the network dynamics is not yet 
sufficiently explored. 



C. External inputs 

Real networks usually have some nodes that respond to 
external inputs. Such an external input to a node can be 
modelled by switching the constant function from 1 to 
or vice versa. The set of nodes that cannot be controlled 
in this way is called the computational core. Networks 
with a higher proportion of C2 functions tend to have 
a larger computational core, since the C2 functions can 
mutually fix or control each other. Investigations of this 
type can be found in [35| . 



D. Evolution of Boolean networks 

Ensembles of networks that are completely different 
from the random ensembles studied in this review can 
be generated by evolving networks using some rule for 
mutations and for the network "fitness" . For instance, by 
selecting for robustness of the attractors under noise, one 
obtains networks with short attractors that have large 
basins of attraction, but that do not necessarily have a 
large frozen core [H, [13, HU ■ 

In another class of evolutionary models, fitness is not 
assigned to the entire network, but links or functions 
are changed if they are associated with nodes that do 
not show the "desired" behavior, for instance if they are 
mostly frozen (or active) , or if they behave most of the 
time like the majority of other nodes [39l. liH, l4l| . 
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E. Beyond the Boolean approximation 

There exist several examples of real networks, where 
the essential dynamical steps can be recovered when us- 
ing simple Boolean dynamics. If a sequence of states 
shall be rcpcatable and stable, and if each state is well 
enough approximated by an "on" -"off" description for 
each node, Boolean dynamics should be a good approxi- 
mation. However, wherever the degree of activity of the 
nodes is important, the Boolean approximation is not 
sufheent. This is the case for functions such as continu- 
ous regulation or stochastic switiching or signal amplifi- 
cation. Clearly, in those cases a modelling is needed that 
works with continuous update functions or rate equations 



based on concentrations of molecules. The different typ es 
of network modelling are reviewed for instance in [421 ] . 
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